#include "../vector.h"
#include "../sys/rand.h"
#include <math.h>
#include <time.h>

static int compare_reals( const void *lhs, const void *rhs )
{
	const Real a = *(Real *)lhs;
	const Real b = *(Real *)rhs;
	return a < b ? -1 : ( b < a ? 1 : 0 );
}

int main( int argc, char *argv[] )
{
	register size_t i;
	mds_srand( time(NULL) );
	/* test 32-bits rand */
	{
		uint32_t *r = NULL;
		size_t    n = 10000;
		r = MDS_CALLOC_AS(uint32_t,n);
		if( r )
		{
			for(  i=0; i < n; ++i )
			{
				r[i] = mds_rand();
			}
			
			MDS_FREE_AS(uint32_t,r,n);
		}
		else {
			mds_perror(stderr, "allocate r");
		}
		
	}
	
	/* test uniform [0:1] rand */
	{
		Real     *r = NULL;
		size_t    n = 20000;
		r = MDS_CALLOC_AS(Real,n);
		if( r )
		{
			Real ave = 0;
			Real sig = 0;
			for( i=0; i < n; ++i )
			{
				r[i] = mds_alea();
			}
			qsort(r, sizeof(Real), n, compare_reals);
			
			for(  i=1; i < n ; ++i )
			{
				ave += r[i];
			}
			ave /= n;
			fprintf( stderr, "uniform ave= %.6f / 0.5\n", ave );
			
			for( i=1; i < n ; ++i )
			{
				const Real d = r[i] - ave;
				sig += d*d;
			}
			sig = sqrt( sig /(n-1) );
			
			fprintf( stderr, "uniform sig= %.6f / %.6f\n", sig, 1.0/sqrt(12.0) );
			
			MDS_FREE_AS(Real,r,n);
		}
		else {
			mds_perror(stderr, "allocate r");
		}
		
	}
	
	/* test normal [0:1] rand */
	{
		Real     *r = NULL;
		size_t    n = 20000;
		r = MDS_CALLOC_AS(Real,n);
		if( r )
		{
			Real ave = 0;
			Real sig = 0;
			for( i=0; i < n; ++i )
			{
				r[i] = mds_normal();
			}
			qsort(r, sizeof(Real), n, compare_reals);
			
			for( i=1; i < n ; ++i )
			{
				ave += r[i];
			}
			ave /= n;
			fprintf( stderr, "normal ave= %.6f / 0\n", ave );
			
			for( i=1; i < n ; ++i )
			{
				const Real d = r[i] - ave;
				sig += d*d;
			}
			sig = sqrt( sig /(n-1) );
			
			fprintf( stderr, "normal sig= %.6f / 1\n", sig );
			
			MDS_FREE_AS(Real,r,n);
		}
		else {
			mds_perror(stderr, "allocate r");
		}
		
	}
	
	
	/* test vectors rand */
	{
		mds_vector *a = NULL;
		mds_vector *b = NULL;
		size_t      n = 10;
		const mds_vector box = { 1,2,3 };
		
		a = MDS_CALLOC_AS(mds_vector,n);
		b = MDS_CALLOC_AS(mds_vector,n);
		
		if( a && b )
		{
			for( i=0; i < n; ++i )
			{
				mds_vector *u = a+i;
				mds_vector *v = b+i;
				mds_vector_rand(v,&box);
				mds_vector_unit(u);
				fprintf( stderr, "u=[%10g %10g %10g], length=%10g\n", u->x, u->y, u->z, mds_vector_length(u) );
				fprintf( stderr, "v=[%10g %10g %10g], length=%10g\n", v->x, v->y, v->z, mds_vector_length(v) );
				
			}
		}
		
		MDS_MEMORY();
		MDS_FREE_AS(mds_vector,a,n);
		MDS_FREE_AS(mds_vector,b,n);
		MDS_MEMORY();
	}
	return 0;
}
